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ABSTRACT 

This paper describes research done to apply the Fractional Step Method to finite-element 
simulations of natural convective flows in pure liquids, permeable media, and in a directionally 
solidified metal alloy casting. 


The Fractional Step Method has been applied commonly to high Reynold’s number flow 
simulations, but is less common for low Reynold’s number flows, such as natural convection in 
liquids and in permeable media. The Fractional Step Method offers increased speed and reduced 
memory requirements by allowing non-coupled solution of the pressure and the velocity 
components. 


The Fractional Step Method has particular benefits for predicting flows in a directionally 
solidified alloy, since other methods presently employed are not very efficient. Previously, t e 
most suitable method for predicting flows in a directionally solidified binary alloy was the 
penalty method. The penalty method requires direct matrix solvers, due to the penalty term. e 
Fractional Step Method allows iterative solution of the finite element stiffness matrices, t ere y 
allowing more efficient solution of the matrices. The Fractional Step Method also lends itse o 
parallel processing, since the velocity component stiffness matrices can be built and solve 
independently of each other. 


The finite-element simulations of a directionally solidified casting are used to predict 
macrosegregation in directionally solidified castings. In particular, the finite-element simulations 
predict the existence of “channels” within the processing mushy zone and subsequently 
“freckles” within the fully processed solid, which are known to result from macrosegregation, or 
what is often referred to as thermo-solutal convection. These freckles cause material property 
non-uniformities in directionally solidified castings; therefore many of these castings are 
scrapped. 


The phenomenon of natural convection in an alloy under-going directional solidification, or 
thermo-solutal convection, will be explained. The development of the momentum and continuity 
equations for natural convection in a fluid, a permeable medium, and in a binary alloy 
undergoing directional solidification will be presented. Finally, results for natural convection in 


a pure liquid, natural convection in a medium with a constant permeability, and for directional 
solidification will be presented. 


INTRODUCTION 

Natural convection occurs in many physical situations. Its most common occurrence arises from 
an unstable thermal gradient in a gas or liquid. Thermally driven natural convection can also 
occur in permeable media. Natural convection can also occur in much more complex situations, 
such as in a liquid metal alloy undergoing directional solidification. Here, the driving force can 
be thermal, but density differences driven by concentration differences in the liquid can also 
occur. 

The ultimate goal of the work described here was to develop a mathematical model and 
subsequently the capability to simulate the fluid flows occurring in a directionally solidifying 
alloy, using the fractional step method. The fundamental equations describing flow in a 
directionally solidifying alloy are very complicated. In addition, implementation of the fractional 
step method is challenging. For this reason, the fractional step method was first applied to the 
equations governing natural convection in a single-phase fluid that can be treated as containing 
only a single constituent. Next, the fractional step method was applied to the equations 
governing natural convection in a permeable medium, also containing a single-phase fluid of one 
constituent. Application of the fractional step method to the simulation of this type of physical 
situation is slightly more complicated than that for a fluid. Finally, the fractional step method 
was applied to the equations governing natural convection in a binary alloy undergoing 
directional solidification. 

This paper is laid out as follows. The section following discusses the phenomena occurring in a 
directional solidifying alloy that contribute to natural convection and the reasons for needing to 
understand the flow. Fundamentals of the fractional step method are discussed next. The 
application of the fractional step method to the equations governing flow in a fluid and in 
permeable media is presented next. Presentation of the fractional step method as applied to the 
equations governing flow in a binary alloy undergoing directional solidification follows. Next, 
results of the simulations for the three different flow situations are presented. Finally, 
conclusions will be discussed. 

DIRECTIONAL SOLIDIFICATION 


DEFECTS IN ALLOYS AND THEIR CAUSES 

Directional solidification is a process by which aircraft engine turbines are produced. Casting 
defects often necessitate that the entire casting be scrapped, which is very costly. A major defect 
that leads to scrapped castings is that of concentration non-uniformities. These casting non- 
uniformities are a direct result of the phenomena occur during the solidification process and 
affect the consistency and integrity of the casting’s structural properties. Two phenomena that 
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produce concentration non-uniformities and affect their severity are microsegregation and 
macrosegregation. Microsegregation is a phenomenon that occurs over a distance approaching 
the size of dendrites in the mushy zone and is dependent on the shape of the dendrites and solute 
diffusion in the liquid and solid [1]. Macrosegregation causes compositional differences to occur 
over a distance approaching that of the size of the casting or ingot, and is normally large 
compared to that caused by diffusion or microsegregation. Consult [1,2,3] for additional 
information. 

During the directional solidification process, the casting is cooled from below. This minimizes 
convection by imposing a stable temperature gradient. Strictly speaking, due to the stability 
impressed on the melt by cooling it from below, thermal convection is suppressed. However, 
depending upon various parameters such as the Solute Rayleigh Number, temperature gradients 
perpendicular to the gravity vector can interact with the solute field to cause sideways diffusive 
instability [4,5]. Even more significant is the phenomenon of solute convection occurring when 
solute rejected at the solidification interface is lighter than the solvent. This combination of 
phenomena is often referred to as thermosolutal (double diffusive) convection. Thermosolutal 
convection can cause severe macrosegregation of material. 

This severe macrosegregation in the liquid phase causes defects to occur. These defects are made 
up of long narrow streaks oriented roughly parallel to the direction of gravity. These streaks, or 
“freckles” have solute concentrations that vary significantly from that of the surroundings [6]. 
These freckles make the casting more prone to fatigue cracking at the interfaces. Therefore, the 
castings are scrapped. Several experimental works with nonmetallic transparent systems [7,8] 
have verified this. 


PRESENT SIMULATION CAPABILITIES 

Previous research, related to the research proposed herein, as described in Reference [9] has 
come a long way in simulation of solidification processes. The models currently used by 
Heinrich et al. [9,10,11,12,13] include interdentritic diffusion, thermosolutal convection, and 
solute and energy conservation. Since the mushy zone is a metallic dendrite assembly of low, 
non-uniform, and anisotropic permeability, it is simulated as a region of anisotropic permeability 
with no predetermined size or shape [14]. 

Presently, the size of the domains that can be simulated is limited because of the large amount of 
computer resources required. This limitation is attributed mainly to two issues: 1) In order to 
adequately capture the diffusion processes occurring in the mushy zone, very fine meshes are 
required, on the order of dendrite spacing, and 2) since permeability theory is used to model the 
varying liquid fraction of the mushy zone, the only numerical algorithms that have been 
successfully employed to solve the momentum equation stiffness matrix are direct solvers. 

The penalty method is presently the algorithm employed most frequently for the solution of the 
solidification momentum equation stiffness matrices. The pressure is not calculated directly in 
the penalty method formulation. Rather, the pressure term is eliminated by equivalencing it to 
the continuity equation (which, numerically, has a value very near zero) times a very large 
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penalty parameter. All velocity components are calculated simultaneously at each time step. 
Therefore, for a three-dimensional finite-element computation domain consisting of N nodes, the 
momentum equation stiffness matrix has a size of order 3N X 3N. The stiffness matrix generated 
by the penalty method requires direct matrix solvers, due to the numerical stiffness introduced by 
the penalty term and due to the permeability term. 

The computer resources required by direct formulations are even more than for the penalty 
method, since the matrices generated for a direct mixed formulation are non-symmetric. The 
pressure is calculated with no approximations; therefore there are four unknowns per node, 
making the stiffness matrix size of order 4N X 4N. Since direct formulations characteristically 
produce very large ill-formed stiffness matrices, direct formulations are normally used only for 
calibration of the penalty method and other methods [9]. 

The Gallerkin Least Squares method was used in the 3D computations described in References 
[13] and [15], and allowed for large domains to be simulated. However, the use of this method is 
restricted and very difficult to use. 

The attractive feature of the fractional step method is that the velocity components and the 
pressure can be de-coupled and solved independently. This leads to stiffness matrices that are 
much smaller and therefore much more manageable. In fact, this allows the solution of flow in 
large domains that were not realistic previously. The fractional step method is described next. 


FUNDAMENTAL PRINCIPLES OF THE FRACTIONAL STEP METHOD 

The fractional step method has been used extensively for solution of the incompressible Navier- 
Stokes equations in Computational Fluid Dynamics (CFD). Chorin [16,17,18] introduced the 
method over thirty years ago and Temam [19] independently came up with similar ideas. A 
general description of its implementation is discussed here. 

Beginning with the dimensionless form of the incompressible Navier-Stokes equation and 
continuity: 


~ V=u.„ = -Vp„, -u. ■ Vu. +f. (1.1) 

At Re 

V-u = 0 0- 2 > 

where the velocity has been discretized in time, the viscous and pressure terms are treated 
implicitly and the convective and forcing term f n are treated explicitly. This equation is split into 
two steps, the intermediate “viscous” step and the projection or “inviscid” step. 
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(1.3) 


U " +1 A = U ® ' VU o +f » 

At Re 

and 

Un+1 ~ U ~— ~-^-V 2 (u n+t -u„+, ) = -Vp n+ , (1.4) 

At Re 

The intermediate velocity u n+1 calculated in equation (1.3) does not necessarily meet continuity 
requirements. This intermediate velocity is corrected in equation (1.4) so that u n+1 , the velocity 
at the end of the time step, does meet continuity. The pressure and subsequently the pressure 
gradient, are calculated by removing the implicit viscous term in equation (1.4), and then taking 
the divergence. Due to continuity, equation (1.2), V -u n+1 must vanish and the pressure is 

obtained from: 


V 2 p n+1 = 


At 


(1.5) 


Removing the implicit viscous term or dealing with it in alternative ways makes the calculation 
of the pressure much more straightforward. The reason most often cited for removing this term 
is that the intermediate velocity profile is projected into an “inviscid” vector space, orthogonal to 
the intermediate “viscous” space [20,21]. Some references treat the viscous term explicitly [22] 
in equation (1.3) so that it does not appear in equation (1.4). However, explicit treatment of the 
viscous term restricts numerical stability, particularly for low Reynolds number flow fields. 
Other references include a more mathematically pure derivation [23,24]. In these references, the 
fractional step method derived above is modified as follows: 

At _2 

The pressure in equation (1.4) is replaced with p - yf — V iff : 


M, +l ~ U n + . 

At 


— -V 2 (u 
Re v 


n+1 


-U 


n+1 





( 1 . 6 ) 


The following terms are removed from equation (1.6) to take on the standard form of the 
projection step: 


U n+I fi »±L = -Vy n4l (1.7) 

At 

At . 

Applying the Laplacian Operator to equation (1.7) and multiplying both sides by — yields an 

equality for the remaining terms in equation (1.6); therefore removing them is mathematically 
correct. 

Equation (1.5) is now a Poisson equation for if/ instead of p: 
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( 1 . 8 ) 


vV„ + , 


V u 


n±l 


At 


and the pressure at the end of the time step is calculated as: 

(L9) 

Whether calculating the pressure from equations ( 1 .6), ( 1 .8), and ( 1 .9) is more accurate than 
calculating it directly from equations (1.3) and (1.5) is application dependent and must be teste 

in practice [23]. 


THE ACCURACY OF THE FRACTIONAL STEP METHOD 

A number of authors have predicted the accuracy of the fractional step method for 
incompressible flow [23,24,25,26]. The general consensus is that the velocity can be determined 
to second-order accuracy. However, there is disagreement on pressure. Some claim the pressure 
can be determined to second order accuracy [27,28], while others claim that the very nature of 
the fractional step method limits the pressure to first-order accuracy no matter what steps are 
taken to improve the method [25,26]. Fortunately, there is consistent agreement on how to 
improve its accuracy. 

The three main methods of ensuring acceptable accuracy for the fractional step method are. 


1) Accurate selection of boundary conditions for the intermediate velocity as suggested by 
Kim and Moin [29]. 

2) Accurate boundary conditions for the pressure by Orszag et al. [30] and E. and Liu 
[21,31]. 

3) Pressure-correction schemes as demonstrated by Van Kan [28], Bell et al. [27], 
Quartapelle [22], and Gresho [32,33]. 


The first two items hint at one of the inherent stumbling blocks of the fractional step method 
[25]. There is only one set of the “real” boundary conditions: the boundary condition for the 
divergence-free velocity. Since the splitting of the Navier-Stokes equations into two steps 
requires two sets of boundary conditions, one must be very careful in setting up these boundary 

conditions. 


The third item substitutes the calculation of the pressure at every time step with an “incremental 
pressure calculation. That is, the pressure is adjusted or corrected at every time step relative to 
the previous time step’s pressure. 


MATHEMATICAL MODEL FOR NATURAL CONVECTION: WITH AND WITHOUT 
PERMEABILITY 
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The development of the equations for the solution of Navier-Stokes equations for natural 
convection in a single-phase, single-constituent fluid and in a permeable medium are very 
similar. The only difference is the permeability term. Therefore, the development of the 
equations for a permeable medium are developed here. The equations for natural convection in a 
fluid can be obtained by removing the permeability term. 


DEVELOPMENT OF DIMENSIONLESS MOMENTUM EQUATION FOR NATURAL 
CONVECTION IN A PERMEABLE MEDIUM 

Beginning with the dimensional momentum and continuity equations for a Newtonian Fluid with 
Constant Density and Constant Viscosity, as found in Appendix C.5 of Panton [34], and with the 
permeability term added: 


P 0 


du' 

dt' 


+u'V'u' 


\ 


= - V'p' + pV'V - n (K')~‘ u' + p i 


) 


( 2 . 1 ) 


V'u' = 0 


where the permeability vector K' is defined from the Darcy equation [35]: 

u' = -— (V'p'-pg) (22) 

P 

The primes indicate dimensional terms. Besides the addition of the permeability term, the 
original Panton equation has been slightly revised ( p 0 on left hand side) to allow using the 
Boussinesq approximation [36]. Equation (2.1) shows the permeability term as a vector to 
account for its anisotropy. For the analysis described in this section, K is assumed the same in 
both the x and y directions, therefore, the permeability term will be shown as a scalar term. 

For convenience, a component of the hydrostatic pressure (based on the reference density) is 
separated from the pressure: 


p' = p'* + Po (g'x x ' + gy / ) (23) 

where p'* is a modified pressure, and the reference density p 0 is defined as the density of the 
fluid at T = T 0 . The gradient of the pressure can be written as: 

V'p' = V'p'* + p 0 g (2-4) 

Insert equation (2.4) into equation (2.1) and divide all terms by p 0 : 
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Now use the Boussinesq approximation for p, p = p 0 [l + fi T (T' -T 0 ')] . The dimensional 
momentum equation (after canceling the g terms) becomes: 


9U / r-i/ / 1 Xl f * 

— — + u V u = V ] 

dt' Po 

Non-dimensionalize according to the following reference variables 


- + u' • V'u' = — — V'p" + v V' V - vK" 1 u' + j8 T (T' - ) g' 

dt' Po 


Characteristic Length H 


Characteristic Velocity V = JPrgtfH 
x y u 

x =h' y= T U= V 


* p I • PoA 

p z=.-L— where P 0 =—^ 
D„ H 


P 0 

V = V7 i Da = Jr 

r r-r; 

AT 


g = 


g_ 

g 


f' . H 

t = — where T = 

t 


Substituting these terms into Equation (2.6): 


. Vu = - 51 Vp + Vhl - Da 1 u + ft g ATT g 


vV 


( 2 . 6 ) 


(2.7) 


H 2 at ’ H H 3 r H z H l 

Substitute the characteristic velocity V into Equation (2.8): 


( 2 . 8 ) 


+ 


L£ + ft gAT „.Vu = -£vp- 
v V /3 T g ATH V 2 .1 + /? ATTg g 

tt2 tt2 


H 


H 


(2.9) 


Divide by /3 T gAT : 


( 2 . 10 ) 


VHD X du 

^ T gAT H 2 at 


+ u-Vu = 




ftgATH 


r v P* 


-V 2 u- , V>/ ^ — Da~'u+Tg 


^ft-gATH 2 ^ T gATH 2 

For convenience, drop the * from p and organize according to the following dimensionless 
parameters: 


p T g NTH 3 

The Rayleigh Number : Raj 


vD t 


( 2 . 11 ) 


77ie Prandtl Number : Pr = — 

LAr 


This results in the dimensionless momentum equation for natural convection in a permeable 
medium: 


1 du „ 

. — + u-Vu = - 

7PrRa T dt PrRa 


L_V P+ P*-V 2 u- Da'u + Tg (2.12) 
Ra T ^ Ra T ^ Ra T 


The dimensionless continuity equation has the same form as the dimensional version: 

Vu = 0 (2-13) 


APPLICATION OF THE FRACTIONAL STEP METHOD TO THE MOMENTUM 
EQUATION FOR NATURAL CONVECTION IN A PERMEABLE MEDIUM 

In order to build stable stiffness matrices for the finite element method, equation (2.12) is 
rearranged to treat the viscous term and permeability term implicitly and all other terms are 
moved to the right hand side to be treated explicitly. The fractional step method requires that the 
transient term be discretized, therefore the subscripts n and n+1 are added, where n refers to the 
value at the previous time step (known) and n+1 refers to the value at the present time step 
(unknown). 


u P +i ~u. 


At^PrRa 1 


Pr 

Ra x 


Da -1 u n+1 - 


f Pr 
Ra T 


* u n+1 = - 


PrRa T 


Vp n+1 -u n • Vu n +T n g (2.14) 


This is rearranged to the following: 
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, 1 — (i + AtPrDa" 1 -AtPrV 2 )u n+1 =— — ^ — v P. + i + f n 


(2.15) 


where the forcing term f„ = -u„ • Vu„ +T„ g + 


u„ 


At Raj 


Use the equality p n+1 = p n +(p n+1 - p n ) to allow calculation of the incremental pressure at every 

time step. Equation (2.15) then becomes: 

1 , _ _ . . . r Vp„ + V (p„ +1 -A,)] 

— (7 + AtPr Da~' -flArPrV )u„ +1 =-^ — 

AtjPrROr V Pr/fa r 


+ f_ 


(2.16) 


Split Equation (2.16) into a intermediate (viscous) step and a projection (inviscid) step, 
respectively: 


====(/ + At Pr Da" 1 -At PrV 2 )u n+1 = 

At yjPr Rdj 


Pr Raj 


Vp n +f„ 


(2.17) 


fL=(l + At Pr Da"' - At Pr V 2 )(u n+1 -u B+1 ) = V (p n+1 - P n ) (2-18) 

At^PrRaj. v pr/te r 

where u , is the solution to the intermediate step and does not meet continuity. Now, take the 

/t+1 

divergence of the projection step, equation (2.18): 


At 


-jL= -V {(/ + At Pr Da" 1 - At Pr V 2 )(u n+1 -u B+1 )} = (p*+ -i _ P * ) (2 ‘ 19) 


Note that at the end of the time step, the velocity must meet continuity, therefore the divergence 
of both u n+1 and V 2 u n+1 are 0: 


^- Pn )--^(l^Da 1 -AtPrV 2 )Vu n+1 (2.20) 

and the viscous term is ignored: 

V 2 (p fl+ ,-pJ = ^^-(^ + AtPrDa" , )v u n+1 (2.21) 

To find the divergence-free velocity at the end of the time step, solve for each component of 
velocity from equation (2.18), again with the viscous term ignored. 
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1 


( 2 . 22 ) 


At 


“„ + i = M « + r 


(l + At Pr Da ~ 1 ) yJPrRa^ dx 


-V ( Pn+1 Pn ) 


V '“ = ’« " (l + A ,La-')>^| (P "'' ‘ P " ] <2 ' 23) 

For a pure liquid, the Darcy term K and its dimensionless counter-part Da approach infinity, 
therefore K~ l and Da approach zero and can be eliminated from the above equations. In order 
to solve the above equations numerically, the finite element method is applied. The particular 
type of finite element method applied is the Gallerkin Method. Consult Heinrich and Pepper 
[37,38]for details on how the finite element method was applied. 


THE FRACTIONAL STEP METHOD APPLIED TO THE MOMENTUM EQUATION 
FOR DIRECTIONAL SOLIDIFICATION IN A BINARY ALLOY 

The above formulation demonstrates the process of applying the fractional step method to a 
single-constituent, single-phase fluid, with and without permeability. Directional solidification 
of a binary alloy is much more complicated. This process includes two constituents, two phases 
and also has anisotropic permeability properties. Developing the formulation here would require 
considerable space. Therefore, only the results of the formulation are presented. Please consult 
[9] for more details. 


The simulation of directional solidification in a binary alloy includes conservation of energy, 
species, mass, and momentum balance. Since the fractional step method affects only the 
conservation of mass and momentum balance, only those equations are presented here. Please 
consult [9,39,40] for more details regarding the equations for conservation of energy and species. 


The dimensionless equation for momentum balance in a directionally solidifying binary alloy is 
[9]: 


+ A Da-'u - _L V 2 u = -0 Vp - -u • Vu 

3t Re Re 0 


B 30 f$ d(f> . 
±L_r u +- t — v-^-+^ 

<p at 3 Re at 


Ra, 


-T+- 


Ra c 


Re 2 Pr Re 2 Sc 


(s L -0 


g 


(3.1) 


and the continuity equation is: 


•u = /3— . 

p at 


(3.2) 


where is the shrinkage coefficient and is defined as P = 


Ps-Pi 

Pi 
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The momentum equation, equation (3.1) is split up into two steps in order to be solved using the 
fractional step method. The dimensionless equation for the intermediate step is: 


— 1 + 
At 


^ At 
Re 


Da 1 


1 


U „ + I =f n -0Vp n 


(3.3) 


where the (explicit) forcing term f n is defined as: 




--u -Vu - 


Mf^+JLv^+0 

<j> at 3 Re at 



(3.4) 


The dimensionless equation for the projection step is: 


_1_ 

At 


Re 


-a..,)=^v(p.„ -p.) (3.5) 


The implicit viscous term is ignored and continuity is imposed by taking the divergence of 
equation (3.5): 


v[oV(p„, 



V U 


n+t 


At 


At at 


(3.6) 


where the dimensionless term a is a function of the fraction liquid <p and the anisotropic Darcy 
term and is a vector quantity: 


o = 




where the individual components of o are as defined by Sabau, et al. [41]: 


(3.7) 


a, = 


1 + 


± 

ift At 

Re Da, 


and 


CT y=- 


1 + 


0At 
Re Da, 


(3.8) 


The dimensionless scaling parameters are: 
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The Reynold's Number: 


Re = 


VH 


The Thermal Rayleigh Number: 
The Solute Rayleigh Number: 
The Prandtl Number: 

The Schmidt Number: 


Ra x 

Ra s 


ftgOH 4 

VD T 

ftgS 0 H 3 

vD s 




(3.9) 


As was done for the case of natural convection in a liquid or a permeable medium, the finite 
element method was applied in order to solve the above equations numerically. 


RESULTS 

Steady-state results are presented for natural convection in a fluid and in a permeable medium. 
Transient results are presented for a binary metal alloy under-going directional solidification. 


RESULTS OF SIMULATING THE MOMENTUM EQUATION FOR NATURAL 
CONVECTION IN A PERMEABLE MEDIUM 

Figures 1-6 show steady-state simulation results of the fractional step method applied to the 
momentum equation for natural convection in a fluid and in a permeable medium. Note that 
Figures 1 — 5 are contour plots of the temperature with the velocity vectors super-imposed. 
Figure 6 is a contour plot of the pressure for the same case as Figure 5. There are 41 nodes in 
each direction in all cases presented. 

The results are presented in dimensionless form, with a dimensionless area of 1 X 1. The 
dimensionless temperature is 1.0 at the bottom and 0.0 at the top. The gravity vector is oriented 
in the minus y direction in all simulations, to create an unstable thermal condition. The gravity 
vector has a magnitude of 9.81 m/s in all simulations. Table I summarizes these results, listing 
the dimensionless parameters defined in equation (2.1 1). Fluid properties are given, along with 
appropriate scaling parameters that correspond to the Rayleigh Numbers and Prandtl Numbers 
listed. In two of the cases listed, fluid properties of air and water at 300 K were used. Consult 
equation (2.7) to determine how reference variables relate to physical dimensions, boundary 
conditions, and fluid properties. 


RESULTS OF SIMULATING THE MOMENTUM EQUATION IN A BINARY ALLOY 
UNDERGOING DIRECTIONAL SOLIDIFICATION 
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Figure 7 shows the rectangular simulation domain of a 0.008 m by 0.03 m rectangle for an alloy 
comprised of 94.2 % Lead, 5.8 % Antimony by weight (Pb-5.8% wt. Sb). A mesh spacing of 
0.0002 m is used, which is on the order of the primary dendrite arm spacings. The initial 
condition consists of a fixed temperature at the bottom boundary (Y’ = 0.0) set to the alloy 
liquidus temperature and an imposed thermal gradient of G 0 = 12,000 K/m (positive with respect 
to Y’) throughout the liquid melt. The liquid melt is of uniform composition S 0 (5.8%). The 
simulated boundary conditions are an imposed thermal gradient of 12,000 K/m at the top 
boundary, and 0.36 K/s cooling at the bottom (the temperature at the bottom boundary is changed 

in time). A uniform body force of 9.81 m/s is imposed in the minus Y’ direction. Figure 7(a) 
shows the three zones that are normally present during solidification: liquid, solid, and mush. 

Figures 8 - 1 1 are simulation results of the fractional step method applied to the momentum 
equation for natural convection in a binary alloy undergoing directional solidification. In this 
case, the boundary condition is continually changing and there is not a steady-state condition. 

The simulated time is 400 seconds in Figure 8-12. 

Figure 8 shows the temperature and pressure, clearly indicating the imposed temperature gradient 
and the resulting pressure gradient. Figure 9 shows a plot of the fraction liquid and the 
concentration of the alloy throughout. The three regions are labeled on the left hand side of 
Figure 9: “Liquid, “Mush” or mushy zone, and “Solid”. At this solidification time, the mushy 
zone is located approximately between Y = 0.007 and 0.010 meters. The “freckles on the 
concentration plot are clearly seen and are a result of thermosolutal convection. The distortion of 
the fraction liquid from a purely vertical gradient is also obvious and is a direct result of the 
thermosolutal convection. Figure 10 shows the same plot of fraction liquid and concentration as 
in Figure 9, with the velocity vectors super-imposed. It is seen that the velocities emanate from 
the regions of highest deviation from the base concentration of 5.8%. Figure 1 1 shows a close- 
up of the velocities and concentration occurring just above the mushy zone. Finally, Figure 12 
shows the velocities and concentration occurring in the mushy zone. The code is achieving 
acceptable continuity when simulating velocities even in the mushy zone, where the magnitude 
of the velocity is several orders of magnitude smaller than the velocities occurring in the all- 
liquid region. 


CONCLUSIONS 

Natural convection is a phenomena common to fluids, permeable media, and also occurs in more 
complex situations such as in a binary alloy under-going directional solidification. Thermo- 
solutal convection in castings undergoing directional solidification leads to concentration non- 
uniformities that necessitate scrapping the castings. The fractional step method has been 
successfully applied to simulation of natural convection in pure fluids, permeable media, and in a 
binary alloy undergoing directional solidification. Its use in simulations of directional 
solidification allows larger domains to be simulated. The results presented here are for two- 
dimensional domains only. Presently, efforts are underway to extend simulation of directional 
solidification with the fractional step method to three dimensions. 
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Table I 

Scaling Parameters, Reference Variables, and Fluid Properties for Figures 1- 6 


Figure 

Ra 

Pr 

Da 

a 

AT 

(K) 

Fluid 

Properties 

H* 

(m) 

v‘ 

(m/s) 

1 

io 5 

1.00 

N/A 

100 

See Note b 

1.0 

0.32 

2 

10 5 

7.00 

N/A 

100 

See Note c 

1.0 

0.26 

3 

10 6 

0.71 

N/A 

85 

Air @ 300 K 

0.05 

0.38 

4 

10 6 

1.0 

io 4 

100 

See Note d 

1.0 

0.31 

5,6 

10 8 

5.83 

io 6 

20 

Water @ 300 K 

0.06 

0.06 


a - ReferenceVariables given in equation (2.7). 

b - The thermal coefficient of expansion is assumed to be —0.0001 K , the kinematic viscosity 

-3 2 . . -3 2 

is assumed to be 10 m /s, and the thermal diffusivity is assumed to be 10 m /s. 
c - The thermal coefficient of expansion is assumed to be -0.0001 K , the kinematic viscosity is 
assumed to be 7 X 10 m /s, and the thermal diffusivity is assumed to be 10 m / s. 
d - The thermal coefficient of expansion is assumed to be -0.001 K , the kinematic viscosity is 
assumed to be 10 3 mis, and the thermal diffusivity is assumed to be 10 m Is . 
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Figure 1: Natural Convection in a fluid, Ra = 10 , Pr — 1. 
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Figure 2: Natural Convection in a fluid, Ra = 10 , Pr — 
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Figure 4: Natural Convection in a Permeable Medium, Ra = 10 , Pr = 1, Da = 10 . 
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Figure 5: Natural Convection in a Permeable Medium with 
Water as the Base Fluid. Ra = 10 , Pr = 1, Da = 10 . 



X 

Figure 6: Natural Convection in a Permeable Medium with Water 
as the Base Fluid. Ra = 10 , Pr = 1, Da = 10 . 
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Figure 7. (a) Basic domain and coordinate system for solidification of a binary alloy, (b) 
Dimensions (m) and mesh spacing of the domain used in the simulations reported here. 
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Figure 8. Temperature and Pressure Results for Pb-5.8 wt. % Sb. 
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Figure 9. Fraction Liquid and Concentration Results for Pb-5.8 wL % Sb. 
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Figure 10. Fraction Liquid and Concentration Results for Pb-5.8 wt. % Sb. 
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Figure 12. Concentration, Velocity Results 
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NOMENCLATURE 


Scalars 

P 
Pi 
Ps 
Po 
Pl 
P s 

0 

p 

v 

¥ 

Da 

Dr 

D s 

g 

G 

H 

I 

K 


shrinkage coefficient due to phase change from liquid to solid 
coefficient of thermal expansion, 1/K 
coefficient of solute expansion 

3 

reference density, kg/m 

3 

density of liquid, kg/m 

3 

density of solid, kg/m 
fraction liquid 

2 

dynamic viscosity, N-s/m 

2 , 

kinematic viscosity, m /s 
pressure function, Pa 

permeability, dimensionless 

2 

thermal diffusivity, m /s 

2 

solute diffusivity, m /s 

2 

value of gravity vector, m/s 
reference thermal gradient, K/m 
reference length, m 
Identity Matrix 

2 

permeability , m 
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p’ pressure. Pa 

p'* modified pressure. Pa 

Pr Prandtl Number 

Re Reynold’s Number 

Ra^, thermal Rayleigh Number 

Ra s solute Rayleigh Number 

Sc Schmidt Number 

S 0 reference solute concentration, % 

t’ time, s 

A t time step size, s 

T’ temperature, K 

AT reference temperature difference, K 

u’ velocity in x direction, m/s 

v’ velocity in y direction, m/s 

V reference velocity, m/s 

Vectors 

g gravity vector, m/s 

n direction normal to boundary 

t’ direction tangent to boundary 

u’ velocity vector, m/s 


